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ABSTRACT 

We study the luminosity function and the radial distribution of satellite galaxies 
within Milky Way sized haloes as predicted in Cold Dark Matter based models of 
galaxy formation, making use of numerical N-body techniques as well as three dif- 
ferent semi-analytic model (SAMs) galaxy formation codes. We extract merger trees 
from very high-resolution dissipationless simulations of four Galaxy-sized DM haloes, 
and use these as common input for the semi-analytic models. We present a detailed 
comparison of our predictions with the observational data recently obtained on the 
Milky Way satellite luminosity function (LF). We find that semi-analytic models with 
rather standard astrophysical ingredients are able to reproduce the observed luminos- 
ity function over six orders of magnitude in luminosity, down to magnitudes as faint 
as My = — 2. We also perform a comparison with the actual observed number of satel- 
lites as a function of luminosity, by applying the selection criteria of the SDSS survey 
to our simulations instead of correcting the observations for incompleteness. Using 
this approach we again find good agreement for both the luminosity and radial distri- 
butions of MW satellites. We investigate which physical processes in our models are 
responsible for shaping the predicted satellite LF, and find that tidal destruction, sup- 
pression of gas infall by a photo-ionizing background, and supernova feedback all make 
important contributions. We conclude that the number and luminosity of Milky Way 
satellites can be naturally accounted for within the (A) Cold Dark Matter paradigm, 
and this should no longer be considered a problem. 

Key words: galaxies: haloes - cosmology:theory, dark matter, gravitation - methods: 
numerical, N-body simulation 



1 INTRODUCTION 

The Milky Way environment provides an excellent labora- 
tory for astrophysics. It has been used extensively in the 
past decades to test theoretical models of galaxy formation. 
In particular, the number density of satellites around our 
Galaxy has long been considered one of the major problems 
for the otherwise quite successful ACDM paradigm. 

About a decade ago, N-body simulations attained suf- 
ficient dynamic range to reveal that in CDM models, all 
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haloes should contain a large number of embedded sub- 
haloes that survive the collapse and virialization of the par- 
ent structure (Klypin et al. 1999; Moore et al. 1999 and more 
recently Diemand et al. 2007). Although the predicted num- 
ber of substructures was in reasonable agreement with ob- 
served luminosity functions in cluster sized haloes, in Milky 
Way sized haloes the number of predicted sub-haloes ex- 
ceeded the number of observed satellites by at least an order 
of magnitude: the known satellite population at that time 
consisted of about 40 satellites with V c > 20 km/s in the Lo- 
cal Group (e.g. Mateo 1998) , while the simulations predicted 
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about 300 sub-haloes with V c > 20 km/s (Klypin et al. 1999; 
Moore et al. 1999). 

Several astrophysical solutions to this problem have 
been proposed. Many authors have pointed out that accre- 
tion of gas into low-mass haloes and subsequent star for- 
mation is inefficient in the presence of a strong photoion- 
izing background, as this background radiation raises the 
entropy of the gas, preventing it from accreting onto small 
dark matter haloes and lengthening the cooling time of that 
gas which has accreted (e.g. Babul & Rees 1992; Efstathiou 
1992; Thoul & Weinberg 1996; Quinn, Katz, & Efstathiou 
1996). Several studies showed quantitatively that this sup- 
pression of gas infall by cosmic reionization could plausi- 
bly reconcile the observed and predicted numbers of "classi- 
cal" Local Group satellites (Bullock et al. 2000, Somerville 
2002, Ricotti, Gnedin & Shull 2002, Benson et al. 2002, 
Read et al. 2006). It was also pointed out that tidal strip- 
ping and heating as satellites orbited in the potential of the 
larger galaxy could cause dramatic mass loss, even decreas- 
ing the circular velocity in the inner parts of the sub-halo 
(Kravtsov, Gnedin & Klypin 2004; Taylor & Babul 2004; 
Zentner et al. 2005). Thus many Local Group satellites may 
inhabit dark matter (sub-)haloes that are much less massive 
than they were at the time that they were accreted by their 
host halo. 

In recent years the Sloan Digital Sky Survey (SDSS: 
Adelman-McCarthy et al. 2008) has changed our view of the 
Milky Way and its environment. The SDSS has made it pos- 
sible to carry out a systematic survey for satellite galaxies, 
which are detectable through their resolved stellar popula- 
tions down to extremely low surface brightness. As a result 
the number of known dwarf spheroidals has doubled in the 
recent past (e.g Willman et al. 2005, Zucker et al. 2006; Be- 
lokurov et al. 2007; Irwin et al. 2007, Gilmore et al. 2007). 
Spectroscopic surveys subsequently measured the velocity 
dispersions of these systems, and confirmed their galactic na- 
ture (Martin et al. 2007, Simon & Geha 2007). This recently 
discovered population of ultra faint satellites has posed new 
challenges for models of galaxy formation and opened the 
possibility to test the ACDM paradigm at very small mass 
scales (e.g. Strigari et al. 2008, Maccio et al. 2009). 

These new observations have made it possible to probe 
the faint end of the luminosity function of Milky Way satel- 
lites, down to luminosities as faint as 100 Lq . Moreover the 
homogeneous sky coverage of the SDSS enables a robust 
determination of the detection limits for faint satellites. Ko- 
posov et al. (2008) provided the first determination of the 
volume corrected Milky Way satellite luminosity function 
down to these extremely faint limits, by assuming various 
simple radial distribution functions for the satellite popula- 
tion and applying the SDSS detection limits. 

In light of the discovery of the new ultra-faint dwarf 
population and the improvements in the numerical mod- 
elling of galaxy formation, it is now timely to revisit the issue 
of whether the basic properties of satellite galaxies around 
the Milky Way, such as their number density, radial distri- 
bution, and mass-to-light ratios, can be reproduced within 
current cosmological ACDM-based models. It is also inter- 
esting to ask what physical processes might plausibly give 
rise to this population of extremely low-luminosity galaxies. 

In this paper we combine merger trees extracted from 
very high resolution N-body simulations with three differ- 



ent semi-analytic model (SAM) codes. These merger trees 
describe the hierarchical assembly of a Milky Way-like halo, 
while the SAMs are used to predict the relationship between 
the dark matter (sub) haloes and observable galaxy proper- 
ties, allowing us to make a direct and detailed comparison 
with observational data. 

The goal of this work is not only to test whether the ob- 
served properties of Milky Way satellites, including the re- 
cently discovered faint population, can be reproduced within 
the ACDM model, but also to understand how and when 
this extreme population formed. We aim to understand how 
various mechanisms (such as SN feedback, cosmic photo- 
ionization, and tidal stripping) may shape the luminosities 
of galaxies populating low mass dark matter substructures 
orbiting around Milky Way-like galaxies. 

The remainder of this paper is organized as follows. In 
Section[2]we describe the numerical simulations. In Section[3] 
we briefly summarize the SAMs used in our study, highlight- 
ing the differences among the models. Section U contains a 
detailed description of the observational data used in this 
work. In Section [5] we compare the luminosity function, and 
radial distribution of simulated satellites with observational 
data. Finally in Section [6] we present our main conclusions. 



2 SIMULATIONS 

The N-body simulations of this study were obtained using 
pkdgrav, a treecode written by Joachim Stadel and Thomas 
Quinn (Stadel 2001). The initial conditions are generated 
with the GRAFIC2 package (Bertschinger 2001). The start- 
ing redshift Zi is set to the time when the standard deviation 
of the smallest density fluctuations resolved within the sim- 
ulation box reaches 0.2 (the smallest scale resolved within 
the initial conditions is defined as twice the intra-particle 
distance). The cosmological parameters are chosen to be: 
«a=0.732, fi m =0.268, fi 6 =0.044, h = 0.71 and <t 8 = 0.9, 
and are in reasonable agreement with the recent WMAP 
mission results (Komatsu et al. 2009). 

We selected four candidate haloes with a mass similar 
to the mass of our Galaxy (M ~ 10 12 M Q ) from an exist- 
ing low resolution dark matter simulation (300 3 particles 
within 90 Mpc) and re-simulated them at higher resolution. 
Our high resolution haloes all have a quiet merging history 
with no major merger after z = 2, and thus are likely to 
host a disk galaxy at the present time (with the excep- 
tion of G3, which we discuss further below). The standard 
high resolution runs are 12 3 times more resolved in mass 
than the initial simulation: the dark matter particle mass is 
nid — 4.16x 10 5 /i _1 Mq, where each dark matter particle has 
a spline gravitational softening of 355 h~ l pc. Some of the 
main properties of the re-simulated haloes are listed in Table 
[1] One of the haloes, namely G3, has a mass greater than 
the expected mass of the MW and has experienced a major 
merger at z=1.5 so it is likely to host an elliptical galaxy. In 
order to check possible resolution effects (especially in the 
construction of the merger tree) we re-simulated one of the 
haloes (namely Gl) with higher resolution (27 3 times with 
respect to the low resolution), with more than 32 million 
particles within the virial radius (GIhr in Table [TJ, reach- 
ing a dark matter particle mass of m = 3.65 x 10 4 /i -1 M o . 

For the purpose of constructing accurate merger trees 
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Table 1. Dark Matter Halo parameters 



Halo 


Mass 
(1O 12 /i" 1 M ) 


Npa.it 

(io 6 ) 


(kpc/h) 


^circ 

(km/s) 


GO 


0.88 


2.12 


197 


178 


Gl 


1.22 


2.93 


219 


188 


G2 


1.30 


3.12 


250 


203 


G3 


2.63 


5.64 


268 


236 


GIhr 


1.15 


31.5 


211 


181 



for each simulated halo, we analyse 53 output times between 
z = 20 and z = 0. For each snapshot, we look for all the viri- 
alized isolated haloes within the high resolution region us- 
ing a Spherical Overdensity (SO) algorithm. We use a time 
varying virial density contrast determined using the fitting 
formula presented in Mainini et al. (2003) . We include in the 
halo catalogue all the haloes with more than 100 particles 
(see Maccio et al. 2007, 2008 for further details on our halo 
finding algorithm) . Our procedure to construct merger trees 
is described in detail in Maccio, Kang & Moore (2009). We 
used all particles within 1.5 times the virial radius of a given 
"root" halo at z — and then track them back to the pre- 
vious output time. We then make a list of all haloes at that 
earlier output time containing marked particles, recording 
the number of marked particles contained in each one. We 
use the two criteria suggested in Wechsler et al. (2002) for 
halo 1 at one output time to be labeled a "progenitor" of 
halo 2 at the subsequent output time. In our language, halo 

2 will then be labeled as a "descendant" of halo 1 if i) more 
than 50% of the particles in halo 1 end up in halo 2 or if ii) 
more than 75% of halo 1 particles that end up in any halo at 
time step 2 end up in halo 2 (this second criterion is mainly 
relevant during major mergers). Thus a halo can have only 
one descendant but there is no limit to the number of pro- 
genitors. On average there are 20,000 progenitors for haloes 
G0-G3, while the number of progenitors for the GIhr run 
is close to 100,000. 

In order to identify subhalos in our simulation we 
have run the MPI+OpenMP hybrid halo finder AHF 
(AMIGA halo finder, to be downloaded freely from 
http:/ /www. popia.ft.uam.es/AMIGA ) described in detail in 
Knollmann & Knebe (2009). AHF locates local overdensi- 
ties in an adaptively smoothed density field as prospective 
halo centres. The local potential minima are computed for 
each of these density peaks and the gravitationally bound 
particles are determined. Only peaks with at least 50 bound 
particles are considered as haloes and retained for further 
analysis. As subhaloes are embedded within their respective 
host halo, their own density profile usually shows a charac- 
teristic upturn at a radius r t r v i r , where r v i r would be 
their actual (virial) radius if they were found in isolation. 
We use this "truncation radius" rt as the outer edge of the 
subhalos and hence (sub-)halo properties (i.e. mass) are cal- 
culated using the gravitationally bound particles inside rt- 

3 SEMI-ANALYTIC MODELS 

We make use of three different semi-analytic model (SAM) 
codes in order to predict the observable properties of galax- 



ies that inhabit the dark matter haloes and sub-haloes iden- 
tified in the TV-body simulations described above (see Baugh 
2006 for a recent review on the semi-analytic approach) . We 
will consider predictions from the most recent implemen- 
tations of three different SAMs, developed independently 
by different groups: (i) the Kang et al. (2005) model that 
has been recently updated in Kang (2008, K08 hereafter); 

(ii) the fiducial model of Somerville et al. (2008, S08 here- 
after), which builds on the original formulation presented in 
Somerville & Primack (1999) and Somerville et al. (2001); 

(iii) MORGANA, first presented in Monaco, Fontanot & Taf- 
foni (2007) and then updated in Lo Faro et al. (2009). Since 
all SAMs assume that DM haloes are the sites where galaxy 
formation takes place and they need a proper description 
of their assembly history, we will use the four merger trees 
extracted from N-body simulations of the G0-G3 haloes (see 
section [2} as a common input. In order to increase the sta- 
tistical robustness of our results, in sec. I5.2l we also consider 
a larger set of realizations of merger trees obtained using the 
extended Press-Schechter (EPS) formalism (e.g. Somerville 
& Kolatt 1999, Parkinson et al. 2008) for K08 and S08 and 
the lagrangian code PINOCCHIO (Monaco et al. 2002) for 
MORGANA. 

All the SAMs considered in this work parametrize in 
different ways the main physical processes acting on the 
baryonic component, such as atomic cooling, cosmic reion- 
ization, star formation, supernovae feedback, metal produc- 
tion and dust attenuation. For sake of simplicity, we will 
discuss here only those processes relevant in shaping the LF 
of MW satellites. We refer the reader to the original works 
for a more detailed discussion on the modeling of physical 
processes (see also Fontanot et al. 2009, for a comparison 
between different SAMs). 

Although our simulations resolve subhaloes we do not 
record the fate of subhaloes in our merger trees or make use 
of this information in the SAMs. When a subhalo is accreted, 
its position is initially either set equal to the virial radius 
of the parent halo at that time (K08 and S08), or extracted 
from a suitable distribution of radial distances (morgana). 
Moreover the orbital parameters (velocity and orbit eccen- 
tricity) for each infalling satellite are randomly selected from 
a distribution motivated by the statistics of satellite orbits 
in cosmological simulations. The dynamical evolution (and 
so the survival probability) of each subhalo is then computed 
by estimating the time required for the subhalo to lose all 
of its orbital energy due to dynamical friction against the 
background DM potential (using updated variants of the 
classical Chandrasekhar formula). 

Each of the models that we have considered applies a 
different set of criteria to determine when satellites are de- 
stroyed by tidal stripping. In the K08 models, a subhalo is 
considered to be tidally destroyed if it either loses more than 
98% of its mass (e.g. Penarrubia et al. 2008) or if its mass 
falls below 6.5 x lO 6 A/0, which is the minimum mass ob- 
served for Milky Way satellites (Strigari et al. 2008). In the 
S08 model, satellites are considered to be tidally destroyed 



1 From this point on, we refer to the DM haloes living within the 
virial radius of larger haloes as "substructure" or "subhaloes", 
while we refer to the all the galaxies except the central galaxy of 
the larger halo as "satellites". 
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when their stripped mass drops below the mass contained 
within a fixed fraction of the halo's original NFW (Navarro, 
Frenk and White 1997) scale radius r s (following Zentner & 
Bullock 2003 and Taylor & Babul 2004). In MORGANA, the 
tidal radius is computed at the first periastron of the satel- 
lite orbit by computing the radius at which the density of 
the unperturbed satellite is equal to the density of the main 
DM halo at the periastron. All the mass (whether dark, stel- 
lar, or gaseous) external to the tidal radius (i.e. at a lower 
density) is then considered unbound. The MORGANA esti- 
mates of the radii of the bulge and disk components, plus 
assumed density profiles for the stars and gas, are used to 
estimate the fraction of the baryonic mass that lies outside 
the tidal radius. 

In all three models the effect of reionization is expressed 
in terms of a "filtering mass" (e.g. Gnedin 2000) . This filter- 
ing mass corresponds to the mass at which haloes will only 
be able to accrete half of the universal baryonic content. 
The fraction of baryons that can be accreted as hot gas is 
parameterized using the following expression (Gnedin 200): 

/b, a cc(2,M vir ) = [1+a26Mp(z)/Mvir ]3 ' W 

where /{, is the universal baryon fraction and M v i r is the 
halo virial mass. The filtering mass as a function of redshift 
Mp(z) depends on the reionization history of the Universe, 
and is parameterized using the fitting formulae provided by 
Kravtsov et al. (2004, but see section [5. 4. II) . 

Massive stars and supernovae may impart thermal and 
kinetic energy to the cold interstellar medium: in the K08 
and S08 models, the rate of reheating of cold gas due to 
supernova feedback is given by an expression of the form: 

SN ( Vdisk \ arh . /rA 
m rh = e y—^J m * ( 2 ) 

where Sq N and a r h are free parameters and m* is the star 
formation rate. K08 and S08 adopt similar values of eg N 
and a r h ~ —2, chosen to reproduce the faint end slope of 
the observed z — global galaxy luminosity function or the 
low-mass end of the stellar mass function. In contrast, the 
MORGANA model adopts a recipe based on the notion of 
a self-regulated feedback loop between star formation and 
supernovae (Monaco 2004), which roughly corresponds to 
Eq = 1 and a r h = in terms of eqn. [2] However, we find 
that in order to reproduce the MW satellite LF with MOR- 
GANA we need to introduce a strong dependence of the mass 
loading factor rj = m r h/m» on the galaxy circular velocitjQ 
(a r h = —4). This suggests that the phenomenological scaling 
in eqn. [2] is a key to the success of galaxy formation models 
in the ACDM context in predicting the MW satellite LF. 

It is worth noting that the ingredients of these SAMs 
have all been developed with much larger galaxies in mind, 
and the models have previously been calibrated mainly 
against observations of relatively luminous galaxies (M* > 
1O 9 M0 or Mv ;$ —16). It is quite unclear whether the stan- 
dard semi-analytic empirical recipes for e.g. star formation 
or supernova feedback should apply in galaxies as tiny as the 
ultra-faint Milky Way dwarfs, which may form out of just 

2 In order not to spoil the agreement of this model for M* > 
IO^Mq galaxies, we retain the original recipe V c > 100 km/s DM 
halos and adopt Eqn. [2]in galaxies with V c < 100 km/s. 



a few molecular clouds. Therefore, it is quite an interest- 
ing experiment to see how well these models perform when 
extended to these very different mass scales. 



4 OBSERVATIONAL DATA 

We test our MW models against observations by focusing 
on two key aspects of MW satellite galaxy properties: their 
luminosity and radial distributions. 

For the luminosity function we use the results of Ko- 
posov et al. (2008; SK08 hereafter). SK08 recently pre- 
sented a quantitative search methodology for Milky Way 
satellites in the SDSS DR5 data and used this method to 
compute detection efficiency maps, which ultimately allowed 
the construction of the first completeness-corrected satel- 
lite galaxy luminosity function (see also Walsh, Willman & 
Jergen 2009). SK08 measured the luminosity function (LF) 
down to My ~ —2, and found that it can be described by a 
single power law of the form dN/dM v = 10 x 10° ' 1(Mv +5) . At 
the very faint end (Mv > —5), in order to compute the com- 
pleteness correction, a radial satellite distribution around 
the host must be assumed; in all the luminosity function 
plots presented in this paper the upper data points (always 
shown as open circles with no error bars) are obtained as- 
suming an isothermal density distribution while the lower 
points (solid circles with error bars) are obtained assuming 
an NFW (Navarro, Frenk and White 1997) distribution. 

We also use a reverse approach to addressing the com- 
pleteness issue by performing the comparison in "observa- 
tional space". Instead of assuming a radial distribution for 
the observed galaxies, we apply the detection criteria of the 
SDSS to our simulations (see section 15.31 for more details) 
and compare directly with the raw data from the SDSS. For 
this comparison we construct a "hybrid" data set. For satel- 
lites brighter than Mv = —9, in order to increase the num- 
ber statistics, we gather together satellites from the Milky 
Way and the Andromeda galaxy (data from Mateo 1998 and 
Metz et al. 2007), and we assign a weight of w = 0.5 to each 
satellite, assuming that current surveys are complete down 
to this limit. For fainter satellites we collect data from Mar- 
tin, de Jong & Rix 2008 (MdJR08 hereafter) and, in order 
to account for the fact that the SDSS surveyed only one- 
fifth of the sky, we set w = 5 for these faint galaxies. The 
adopted Mv threshold for splitting the observational sam- 
ple is justified by the low luminosity of all newly discovered 
satellites both around the Milky Way and the Andromeda 
galaxy (McConnachie et al. 2008). 

In addition, using the same data set described above, 
we compute the cumulative radial distribution of satellites 
(i.e. the number of satellites within a given distance from 
the Sun). Distances for bright satellites are taken from Metz 
et al. (2007), while we use results compiled in MdJR08 for 
faint galaxies (Mv > —9). In computing the radial distribu- 
tion we assign to each galaxy the same weight adopted for 
the LF. 

5 RESULTS 

In this section we present results for both the SAMs and nu- 
merical (dissipational and dissipationless) simulations and 
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compare them with the observational data set described 
in Section [4] First we compare the dynamical evolution of 
satellites in SAMs and in N-body simulations, and then we 
present results for the Luminosity Function (LF) of simu- 
lated satellites and analyze the importance of different phys- 
ical processes (e.g. reionization, stellar stripping and super- 
nova feedback) in shaping the LF. We then present results 
for the radial distribution and compare them with observa- 
tions. 



5.1 Dynamical evolution of sub- haloes in SAMs 
and N-body simulations 

In the SAMs investigated here, we chose to make use of N- 
body based merger trees for "isolated" (or distinct) haloes 
only, and to model the dynamical evolution of satellites 
semi-analytically (see Section [3}. Thus, there will not be 
a one-to-one correspondence between the masses or posi- 
tions of sub-haloes at z = in the SAMs and in the actual 
N-body simulations. In this section we check that the sta- 
tistical distributions of subhalo masses and radii predicted 
by the semi-analytic models are in agreement with those of 
subhaloes identified in the N-body simulations. This com- 
parison is only possible for the K08 and S08 models, since 
the MORGANA model does not explicitly follow the dynamics 
of dark matter substructures. Figure[T]shows the cumulative 
subhalo mass function from N-body and SAMs for our four 
Merger tree G0-G3. The GIhr halo is shown for the K08 
model only. The SAM results are obtained by averaging over 
ten different realizations of the random orbit selection pro- 
cess and are truncated at the N-body mass resolution limit. 
A simple way to quantify the agreement between two dis- 
tributions is to perform a Kolmogorov-Smirnov (KS) test 
(Press et al. 1992). In the following we will quote as KS re- 
sults, the probability that two distributions are drawn from 
the same parent population. K08 and S08 SAMs give an 
average KS probability of 88% and 90% respectively, when 
compared to the subhalo mass function from the N-body 
simulations, showing thus good agreement with the numer- 
ical results. 

The radial distribution of satellites is a key piece of in- 
formation for deciding whether a satellite will be detectable 
in the SDSS survey. It is therefore important to also check 
that the SAMs correctly predict the radial distribution of 
subhaloes within the parent halo. Figure [2] shows the radial 
number density of subhaloes (without including the central 
galaxy); while on average there is good agreement between 
SAMs and N-body results for both models (KS test results: 
87% and 88% S08 and K08 respectively), there is a sys- 
tematic off-set between K08 and S08, especially at distances 
< 100 kpc/h. The reason for this off-set resides in the dif- 
ferent fitting formulas for tidal destruction and dynamical 
time implemented in the two models, but is nevertheless too 
small to affect our analysis. 

In the semi-analytic models used here, the initial orbit 
eccentricity of an infalling subhalo is randomly selected from 
a distribution motivated by the results from cosmological N- 
body simulations. We checked that the scatter due to differ- 
ent realizations of the orbit distribution is fairly small, and 
mostly effects the less numerous, massive satellites which are 
not the focus of our study. 
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Figure 1 . Comparison between the subhalo mass function at z — 
(within /? v i r ) in the N-body simulations (solid black line) and the pre- 
dictions from the semi-analytic models of K08 (red dashed line) and S08 
(blue dotted line) . Each panel shows results for a different dark matter 
halo. The (black) thin solid line in the Gl panel shows the results of the 
K08 SAM when applied to the GIhr halo. 
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Figure 2. The cumulative number density of satellites as a function of 
radius from the center of the parent galaxy. The results from the N-body 
simulations are shown by the solid (black) line, the shaded (grey) area 
shows the la scatter over five different realizations of the satellite orbit 
distribution. Results from the K08 and SOS SAMs are shown by the (red) 
dashed line and (blue) dotted line respectively. Each panel shows results 
for a different dark matter halo. 
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5.2 Satellite Luminosity Function 

We now compute satellite LFs using the four merger trees 
obtained from the N-body simulations as common input for 
our semi-analytic models; this allows us to isolate the impact 
of the different physical ingredients in the SAMs because the 
DM halo merger histories are exactly the same in all three 
models. Figure|3]shows the predictions from our three SAMs, 
adopting z r — 7.5 as the redshift of reionization. We plot all 
satellites within R = 280 kpc in order to be consistent with 
the SK08 data set. Solid lines show the mean of satellite 
distribution and the (grey) shaded area shows the la Poisson 
scatter around that mean. 

All of the models predict that the number of satel- 
lites brighter than My = —3 is of the order of ~ 100 and 
could easily be several hundreds (e.g G3). This is in agree- 
ment with recent estimates of the number of observed satel- 
lites obtained by several different approaches (e.g. Tollerud 
et al. 2008; Madau et al. 2008; Koposov et al. 2009). 

All the models show that the total number of satellite 
galaxies depends on the host halo mass. The LF of the GO 
halo (the least massive one, see Table [TJ is almost flat and 
has a lack of satellites at the faint end compared to the MW 
data. On the other hand, in the most massive halo (G3), the 
SAMs predict more satellite galaxies at all luminosities than 
are observed in the MW. This trend between halo mass and 
the LF does not depend on the SAM used to populate dark 
matter substructures and indicates that the total mass of the 
dark matter halo has a quite strong influence in determining 
the normalization (and shape) of the satellite LF. 

All three semi-analytic models considered in this work 
are able to do a reasonably good job of reproducing the 
observational data. The K08 and S08 models (upper and 
middle panels) in figure [3] quite successfully reproduce the 
observed LF for satellite galaxies over the entire luminosity 
range —2 > Mv ^ —16. The MORGANA model (lower panels 
in figure |3J, tends to predict slightly more satellites at inter- 
mediate luminosity ( — 10 > Mv > —12) but it is nonethe- 
less in good agreement with the observations especially if 
an NFW distribution is assumed for observed satellites. In 
the Gl panel for the K08 model we also show the LF ob- 
tained using the GIhr merger tree (dotted black line). We 
see that the higher resolution merger tree produces almost 
indistinguishable results. 

In some cases, certain models and certain haloes show a 
dearth or even absence of the most luminous satellites (e.g. 
S08 and K08 GO, G2 in all models). This should not be a se- 
rious cause for concern, because the variance in the number 
of massive subhaloes hosting these luminous satellites is very 
large and depends on the detailed merger history of the halo. 
Moreover, the predicted number in the SAM is very sensi- 
tive to the random orbit chosen, as discussed in Section [571] 
To assess this issue in a more robust way, we need to in- 
crease the number of merger tree realizations of MW mass 
haloes. We do this using semi-analytic merger trees instead 
of the merger trees extracted from N-body simulations. For 
each model we generated 40 merger trees for DM haloes in 
the estimated mass range of the Milky Way dark matter 
halo, (0.8 - 1.2) xlO 12 /!,~ 1 M (Klypin, Zhao & Somerville 
2002). Each of the SAM codes has its own algorithm for gen- 
crating merger trees: the K08 and S08 model use different 
implementation of the EPS algorithm, while MORGANA uses 



the PINOCCHIO code. Figure [4] shows the averaged luminos- 
ity function for the three semi-analytic models (we tested 
that the sub-halo mass function from the EPS/pinocchio 
trees is in agreement with the one extracted from the Nbody 
simulations) . 

The K08 and S08 model are, again, in good agreement 
with the observational data: they are able to fit the Milky 
Way luminosity function in the range — 15 < My < —2 (KS 
test results: 90.9% for S08 and 90.5% for K08). The K08 
model shows a small deficit of satellites at brighter magni- 
tudes (especially if compared with S08 and MORGANA) but 
this occurs where the number of observed satellites has a 
large error bar due to poor number statistics. Our models 
produce much better agreement with the number of lumi- 
nous satellites with My < — 15 relative to the predictions of 
Benson et al. (2002). Because of the large number of differ- 
ences between the Benson et al. models and those considered 
here, we can only speculate on the source of this difference. 

It is also interesting to note that the K08 and S08 mod- 
els suggest that the LF of ultra-faint satellites has a down- 
ward kink below about Mv ~ —6, in better agreement with 
the SK08 observational results adopting an NFW, rather 
than isothermal, radial density distribution for the satel- 
lites. An NFW distribution for satellites is predicted by hy- 
drodynamic simulations (Maccio et al. 2006). The compari- 
son between the MORGANA model and the observational re- 
sults can only be performed down to Mv = —5 in this case. 
This because the PINOCCHIO code has never been tested on 
such small scales (e.g. Li et al. 2006) and we did not feel 
confident in using merger trees with a mass resolution be- 
low ~ 1O 8 /i~ 1 M . In the tested range for My the mor- 
gana model is also in good agreement with observations 
(KS test: 90.2% ). 

In figure [5] (lower panel) we plot the stellar mass and 
luminosity of galactic satellites versus their dark matter 
subhalo mass; results are shown for the K08 model for 
all haloes, and are similar for the other models. For this 
comparison we used both the present dark matter mass 
(Mdm{z = 0)) and the mass of the subhalo at the time 
of accretion (Mdm{zo,cc))- The difference between the two 
reflects the effects of tidal stripping on the dark matter 
substructure. The correlation between the present day dark 
matter mass and luminosity is quite broad and, for low lu- 
minosities, Mdm{z = 0) at a fixed luminosity spans almost 
3 orders of magnitude. This is because tidal stripping of the 
dark matter subhalo washes out the initial correlation be- 
tween luminosity and MoM{zacc)- The same applies to the 
comparison between the dark and stellar mass of galactic 
satellites, as shown in the upper panel of figure [5] 

Recently, Strigari et al. (2008) pointed out that one of 
the curious properties of the newly discovered population 
of faint satellites is that over four orders of magnitude in 
luminosity, these objects seem to contain a nearly constant 
total mass within a radius of 300 pc. If we focus on the 
faint population in figure [5] with —3 < Mv < —10, we can 
see that the huge scatter in DM mass at fixed luminosity 
or stellar mass provides a partial explanation for this ap- 
parent "common" mass scale for the faint satellites. Maccio, 
Kang & Moore (2009, see also Li et al. 2009) investigated 
this in more detail, and presented a direct comparison be- 
tween the predicted luminosity and the mass within 300 pc 
for faint satellites in our simulations. They argued that the 
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Figure 3. The Milky Way satellite luminosity function predicted by our three semi-analytic models (SOS, K08 and MORGANA from top to bottom) using the 
GO-3 halo merger trees (from left to right), with an assumed rcionization rcdshift z r = 7.5. The median of the satellite distribution is shown by the solid line, 
while the shaded area represents the la Poisson scatter around the mean. Observational data arc taken from SK08 under the assumption of two different radial 
distributions of satellites, NFW-likc (solid circles with error bars) and isothermal (open circles). The arrows on error bars indicate that there is only one galaxy 
in that particular bin, and so the Poisson error is formally 100%. The dotted line shows a single power law fit to the data: dN/dM v = 10 X 10 01 < M V+ 5 ). In 
the Gl panel of the K08 model results for the GIhr run arc also shown as a dashed (black) line. 



inner profiles of haloes that are initially very concentrated 
are less effected by tidal heating than haloes that are less 
concentrated, so that the mass within 300 pc is reduced 
for more massive subhaloes (which are less concentrated) 
relative to less massive ones. When they corrected for this 
concentration-dependent modification of the inner density 
profile, they found that the Strigari et al. (2008) results are 
quantitatively reproduced by our simulations. 



5.3 Luminosity Function in the Observational 
Plane 

The luminosity function of SK08 has been determined under 
certain assumptions for the radial distribution of satellites 
around our Galaxy. It is also interesting to apply the obser- 
vational selection criteria to our simulations and compare 
"in the observational plane", i.e. with the raw data from the 
SDSS without completeness corrections applied. To make 
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Figure 4. Satellite Luminosity function predicted by SAMs using semi- 
analytic merger trees (EPS and PINOCCHIO, see text for details). Symbols have 
the same meaning as in figure^ The upper right panel shows the average 
satellite LF for the three semi-analytic models together. MORGANA results 
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such a comparison we have applied to our sample of satel- 
lites a visibility criterion, in order to determine if a given 
satellite would be detected in the SDSS. We assume that all 
satellites brighter than Mv — — 10 would be visible and in- 
cluded in the SDSS sample. For fainter satellites, we adopt 
a criterion based on both satellite distance and luminosity 
(see SK08 for more details): 



log(i?/kpc) < 1.04 - 0.228 x My 



(3) 



In the above formula the distance R is measured from the 
Sun and not from the center of the galaxy. In order to 
convert our galacto-centric distances into helio-centric dis- 
tances, we assume the Sun to be located at 8 kpc from the 
center of the galaxy (8,0,0); moreover, since for each satellite 
galaxy we only know its distance from the galactic center, 
we randomly assign position angles (4>, 8) to each of them 
(see also Tollerud et al. 2008). We exclude from the com- 
parison galaxies more distant than 280 kpc. Moreover since 
the SDSS covers only approximately one fifth of the sky we 
randomly select 1/5 of our satellites. We then average over 
100 different realizations of this random sampling. Obser- 
vational data for the recently discovered SDSS satellites are 
taken from MdJR08. Figure[6]shows the comparison between 
the observations and the luminosity functions obtained with 
the K08 model (the S08 model gives very similar results and 
this comparison is not possible for the MORGANA model, be- 
cause it does not provide the distance of satellite galaxies 
from the main halo). The direct comparison with the ob- 
servational data confirms our previous results on the lumi- 
nosity function. Our models are able to reproduce the data 
for haloes G0-G2, while halo G3 slightly overproduces the 
number of faint satellites. The agreement between the data 
and models implies that the distance-luminosity relation of 
our satellites is similar to the observed one. 

The application of the selection criteria of the SDSS 
to our simulated data also allows us to compare the radial 
number density of satellites in our models vs observations. 
Results are shown in figure [7] From this figure we see that 
our simulations reproduce both the observed slope and nor- 
malization of the satellite radial distribution for GO and G2 
(KS test results: 96% and 91%) while the agreement is less 
good for Gl and G3 (KS: 76% and 77%). The difference 
for G3 can be ascribed to the overall higher visible satellite 
number (e.g. fig. [3}. One possible explanation for Gl can be 
related to its higher formation redshift (half mass in place) 
than the other two galaxies. This implies that subhaloes will 
have, on average a higher accretion redshift, and thus have 
more time to sink to the center. 

Finally in figure [8] we compare the radial distribution 
of DM substructures and "observable" satellites in the K08 
model (obviously, observability depends on many factors, 
but in this context we simply consider all satellites with 
My < —1 to be "observable"). In the upper panel we show 
the number density radial distribution of all "observable" 
satellites (M v < -1), of faint satellites (-9 < Mv < -1) 
and of classical satellites (Mv < —9). The number- weighted 
distribution of observable satellites (which are dominated by 
the much more numerous faint population) traces the sub- 
halo distribution (but is down by a factor of ~ 2) at small 
radii (within R/R V [ T ~ 0.2), but flattens relative to the sub- 
halo distribution at large radii. This implies that "observ- 
able" sub-haloes are more concentrated near the large galaxy 
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Figure 6. C oraparison between the "raw" SDSS data (uncorrected for 
completeness) and model predictions with the SDSS visibility criteria ap- 
plied to the simulations (see text for details). The observational data arc 
shown by the (black) dashed line, the K08 model is shown by the solid (red) 
line, and the shaded area shows the 1<t Poisson scatter around the mean 
value. 
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Figure 8. Upper panels: number- weighted cumulative fractional radial 
distribution of satellites; the black line shows all substructures (dark and 
observable), the (red) dotted, dashed and dot-dashed lines are for the "ob- 
servable" satellites (My < —1), the faint satellites ( — 9 < My < —1) and 
the classical satellites {{M\r < —9), respectively. Lower panel: same as up- 
per panels but with weighted by satellite mass. The results were obtained 
by averaging over the four haloes G0-G3; only the K08 model is shown. 
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Figure 7. The f raction of satellites as a function of the distance from the 
central galaxy. Observational data arc shown by the (black) dashed line 
with error bars (representing the Poisson noise). The K08 model is shown 
by the solid (red) line. 



than the overall population of subhaloes (see also Kravtsov 
et al. 2004). When satellites are weighted by their mass 
(lower panel of the same figure), the distribution is domi- 
nated by the classical (bright) satellites, and is almost iden- 
tical to the mass-weighted distribution for all sub-haloes. 
Thus, the different radial distribution of "observable" and 
dark satellites is due to the suppression of star formation 
in low-mass haloes due to cosmic reionization and feedback 
processes. 



5.4 Physical Processes that Shape the Satellite 
Luminosity Function 

We have shown that it is indeed possible to reproduce the 
observed Milky Way satellite luminosity function within the 
ACDM model; we now investigate the role of various phys- 
ical processes in shaping the LF in our theoretical models, 
with a focus on the origin of the newly discovered ultra- 
faint satellite population. There are several possible physical 
origins for the ultra-faint satellite population: it can origi- 
nate from (i) object that formed in haloes with T < 10 4 K 
via H2 cooling (e.g. Salvadori & Ferrara 2009); (ii) haloes 
with T > 10 4 K that were inefficient at accreting (hot) gas, 
because of photoionization; (iii) haloes with T > 10 4 K in 
which star formation was inefficient because of strong super- 
nova feedback; (iv) objects that originally had larger stellar 
masses but have experienced significant stellar stripping. 

In the following sections we will investigate scenarios 
(ii)-(iv). The first scenario, in which the ultra- faint dwarfs 
form via H2 cooling in haloes with T < 10 4 K, is not ac- 
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counted for in our models. It is likely however that some 
ultra-faint satellites could form in this way, and this could 
help in explaining the small gap between the theoretical pre- 
dictions and observations for My > —3 (e.g. figure 



5.4-1 Reionization 

We test the effect of our adopted parametrization of cosmic 
reionization on our results, by both varying the reioniza- 
tion redshift within our reference model based on Kravtsov 
et al. (2004) and by applying a simple modification to 
this model to take into account recent results based on 
high resolution hydrodynamical simulations (e.g. Okamoto 
et al. 2008; Hoeft et al. 2006). 

The redshift at which reionization occurs is still quite 
uncertain, but it is bracketed in the range 7 < z r < 15 
(3<r range from Komatsu et al. 2009); moreover, due to the 
fact that reionization proceeds in an inhomogeneous way, 
the actual redshift of reionization for the Local Group could 
substantially differ from the average reionization redshift of 
the Universe (Weinmann et al. 2007). Figure [9] shows the 
impact of varying the reionization redshift on the Gl lumi- 
nosity function for the K08 model using the standard reion- 
ization parametrization (S08 and MORGANA show a sim- 
ilar trend). Without any suppression of gas accretion due 
to reionization, the simulated LF contains too few satellites 
fainter than My = — 5 and too many with My ~ — 9. This 
is because, in absence of reionization, hot gas can cool very 
efficiently via atomic cooling, and every halo can transform 
a large fraction of its gas content into stars before SN feed- 
back shuts star formation down. When the effect of reioniza- 
tion is taken into account, the amount of gas available for 
cooling and star formation is reduced in low mass haloes, 
and many galaxies are shifted from intermediate luminosi- 
ties (—15 < My < —6) to low luminosities (My > —6), 
producing a luminosity function that is close to a power-law 
and is in good agreement with the data. It is also interesting 
to note that the LF is almost insensitive to the redshift of 
reionization (solid line shows results for z r = 7.5, dashed 
line for z r = 17); this is in agreement with earlier results 
obtained by Kravtsov et al. (2004). 

In addition to uncertainty about the redshift of reion- 
ization, there is still a debate about the value of the char- 
acteristic mass, Mf, below which galaxies are strongly af- 
fected by photoionization. Okamoto et al. (2008, see also 
Hoeft et al. 2006), using hydrodynamical simulations, re- 
cently suggested that the actual value of Mf can be signifi- 
cantly lower than values previously obtained (e.g. by Gnedin 
2000). In order to explore the implications of their results, 
we introduce a factor 7 that multiplies the original value of 
Mf(z) as derived by Gnedin (2000, see eq.[Q. We have used 
two constant values for 7, namely 0.2 and 0.5, and a red- 
shift dependent expression 7(2) = (1 + 2) ' /11.8, derived 
from figure Bl of Okamoto et al. (2008). The results for these 
three modified models are shown in figure[l0]for the Gl halo 
and for the K08 model (results from the other models are 
similar). We see that the reduced value of Mf has the effect 
of increasing the number of galaxies with My ~ 10, creating 
a bump in the LF, with respect to the standard case, sim- 
ilar to what we saw in figure [9] for the no reionization run 
(dotted line). This is true especially for strong suppression 
of photo-ionization as in the 7 = 0.2 and 7(2) cases. Nev- 
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Figure 9. Satellite luminosity function of the Gl halo in the K08 model 
for three different reionization rcdshifts. The solid line represents our fidu- 
cial model with z r — 7.5, the dashed line is for z r — 17, and the dotted 
line is for a model with no reionization (but including all other kinds of 
feedback). The shaded area represents the lcr scatter around the mean of 
the z r — 7.5 model. 



ertheless it is still possible to reconcile the simulated LFs 
with the observational data by increasing the reionization 
redshift, as shown in the right panel of figure [lOl where we 
use z r = 11. 

In summary, it is interesting to note the interplay be- 
tween the strength of the suppression of gas accretion due 
to photo-ionization (as reflected in the filtering mass scale) 
and the redshift of reionization. In our standard models, in 
which the suppression is relatively strong (Mf is large), we 
find a weak dependence of the predicted LF on the adopted 
value of the reionization redshift, while in models with a 
lower overall normalization of the filtering mass, we find a 
stronger dependence on the redshift of reionization. 



5-4-2 Stellar Stripping and Tidal Destruction 

In our models, we distinguish between "stellar stripping" 
and "tidal destruction" . In the former, we track the amount 
of stellar material that is stripped from the galaxy as it 
orbits within the parent halo. In the latter, we assume that 
a satellite's baryonic mass is unaffected until its total mass 
is stripped by a critical amount, at which point the satellite 
is simply removed. Obviously, these are aspects of the same 
physical process and it is somewhat artificial to distinguish 
between them. We do so simply to illustrate the sensitivity 
of our results to different implementations of this physical 
process in models. 

The MORGANA model allows for the modelling of stel- 
lar stripping as a satellite galaxy orbits around the parent 
galaxy (see section [3] for more details) , while the K08 and 
S08 models assume that the stellar content of a satellite 
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Figure 10. Satellite luminosity function for the Gl halo in the K08 
model, for different paramctrizations of cosmic reionization. The solid line 
shows results for the standard reionization model (7 — 1, see text for the 
definition of 7) based on Kravtsov et al. (2004). The dotted and short 
dashed lines show results for a filtering mass Mp reduced by 80% and 50% 
respectively. The long dashed line is for a model with a rcdshift depen- 
dent expression for the modification of Alp (see text for more details). The 
shaded area represents the lcr scatter around the mean of the 7 — 1 model. 
The left upper and lower panels arc for z T — 11 and 7.5 respectively. 



remains unchanged until the dark matter halo is stripped 
beyond a certain critical point, at which point the galaxy is 
destroyed completely. 

First we investigate the possible effect of stellar strip- 
ping on the predicted satellite LF using the MORGANA code. 
We compare three different MORGANA runs with no, stan- 
dard (moderate) and high stellar stripping; this latter case 
is obtained by increasing by a factor of three the fraction of 
stripped stellar mass with respect to the standard run (i.e. 
we remove from the satellites three times the mass that is 
beyond the tidal radius) . Results for the luminosity function 
are shown in figure [TT1 In the case of standard stripping the 
average stripped stellar mass is of the order 5% — 10% with 
no mass dependence. Comparing the LF of this case with 
the no-stripping case, it appears that stellar stripping is al- 
most negligible for satellites brighter that My = —5 and 
it marginally affects fainter satellites. When the strength of 
the stellar stripping is increased the mass loss is, of course, 
more important and as much as 40% of the stellar mass can 
be stripped, with a strong dependence on the orbital param- 
eters. However, even in this case the overall effect on the LF 
is relatively small. According to these results, stellar strip- 
ping is not one of the most important processes that shapes 
the satellite LF or produces ultra-faint satellites. 

On the other hand we find that the modelling of tidal 
destruction of satellites does have a significant effect on the 
satellite LF. If we neglect tidal destruction, we find many 
more low-mass subhaloes than are seen in the N-body simu- 
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Figure 11. The effect of stellar stripping on the satellite luminosity func- 
tion. Results from the MORGANA model are shown for different levels of 
stellar stripping (no stripping, standard, and high). The solid line shows 
the moderate (standard) stripping case (as shown in figure^; dotted and 
dashed lines show the no-stripping and strong stripping models respectively. 



lations, and we would also predict many more faint satellites 
than are observed in the MW. Following Zentner & Bullock 
(2003) and Taylor & Babul (2004), the S08 model consid- 
ers a satellite to be tidally destroyed when the mass of the 
dark matter sub-halo has been stripped to a value less than 
or equal to the mass within /dis^s, where r s is the halo's 
original NFW scale radius. Zentner & Bullock (2003) adopt 
/ dis = 1 while Taylor & Babul (2004) adopt / dis = 0.1. Fig- 
ure [12] shows the effect of varying the parameter /dj s in the 
S08 model; we see that this can change the number of faint 
satellites by as much as a factor of ten. We find good agree- 
ment with the observed MW satellite LF for /dis ~ 0.1 — 0.5. 



5.4-3 Supernova Feedback 

Feedback from supernovae is believed to be an important 
mechanism for regulating star formation in low mass galax- 
ies, and it plays a primary role in shaping the LF in semi- 
analytic models. It is then important to disentangle its ef- 
fect from the effect of cosmic reionization discussed in the 
previous section. To investigate this, we first compare our 
reference model with both cosmic reionization and SN feed- 
back, and a run with only cosmic reionization (no SN feed- 
back). The results for the K08 model are shown in figure [HJ1 
the S08 model shows similar behaviour. In the absence of 
stellar feedback, the SAM predicts a deficit of faint satel- 
lites (My > —5) and an overabundance of intermediate 
luminosity satellites ( — 15 < My < —10) when compared 
with the observations. The inclusion of SN-driven outflows 
again moves intermediate luminosity galaxies into the ultra- 
faint regime by removing a significant fraction of the gas and 
thereby suppressing star formation. 
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Figure 12. The effect of varying the tidal destruction parameter /dis 
(see text) on the predicted MW satellite LF in the S08 model, for the Gl 
halo. 
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Figure 13. The effect of supernova feedback on the satellite luminosity 
function. Results from the K08 model with SN feedback switched on (solid 
line with shaded area) and off (dotted line) arc shown. 



As we already discussed in Sec. [3] the original version 
of the MORGANA model, implementing a recipe for super- 
nova feedback derived from Monaco (2004), does not pro- 
duce good agreement with the MW satellite LF. In order 
to reproduce the observed MW satellite LF, we found it 



necessary to modify the stellar feedback modeling using an 
approach similar to K08 and S08. 

We conclude that strongly differential SN feedback (in 
which the outflow rate relative to the star formation rate 
is much higher in low-mass galaxies) plays a key role in 
reproducing the MW satellite LF in these models. 



5.4-4 Satellite Strangulation 

There is another effect that in our models can cause rela- 
tively high-mass haloes (T > 10 4 K) to be inefficient at form- 
ing stars and thus to produce ultra- faint satellites. This ef- 
fect is a result of the fairly standard assumption in SAMs of 
"satellite strangulation" , namely that the hot gas halo that 
could provide new gas to a galaxy is stripped immediately 
when a galaxy becomes a satellite in a larger halo. Thus, 
haloes that are accreted by the parent halo soon after they 
crossed the threshold for atomic cooling (T > 10 4 K) and 
thereafter are starved of any new gas cooling or accretion 
can have very low stellar masses. 

In order to explore the importance of this effect, we 
look for a correlation between satellite luminosity and the 
timescale ts, which we define as: 



fa 



TS 



■ tto 



tc 



(4) 



where the formation time (fform) is defined as the time at 
which the halo reaches a virial temperature T > 10 4 K and 
can first begin to cool, i a cc is the time at which the galaxy is 
accreted by the parent halo and becomes a satellite, and the 
cooling time (f CO oi) is the time needed for the gas to radiate 
away all of its energy via atomic cooling, computed at z = 
Zacc (the standard definition of cooling time used in SAMs). 
Satellites with rg < 1 may have been "strangulated" before 
they were able to cool a significant fraction of their gas. 
Figure [T4l shows that some (about 15 %) of the faint satellites 
have ts < 1 and therefore may have been impacted in part 
by this effect (though photo-ionization and SN feedback may 
still play a role in these objects as well). 



6 DISCUSSION AND CONCLUSIONS 

In the last few years, a new population of ultra-faint dwarf 
satellite galaxies has been discovered around our Galaxy. 
Given these new observational data, it is timely to re- 
visit the long standing problem of the number of satel- 
lites around Milky Way-like dark matter haloes as predicted 
in the ACDM scenario. We address this issue by combin- 
ing high resolution N-body simulations with three differ- 
ent semi-analytic models of galaxy formation. Four high 
resolution N-body simulations are used to create detailed 
merger trees that represent the assembly history of a Galac- 
tic dark matter halo. These merger trees are then used as 
common input for three SAMs for galaxy formation, namely 
the MORGANA model (Monaco et al. 2007), the Somerville 
et al. model, (S08, Somerville et al. 2008) and the Kang et 
al. model (K08, Kang et al. 2008), to study the expected 
abundance and properties of satellite galaxies in the Local 
Group. 

Because the SAMs do not use the explicit information 
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Figure 14. Relation between satellite luminosity and tcj (defined accord- 
ing to eqf^ for the KOS model. Different symbols show different values for 
the total mass of the satellite at the time of accretion, as indicated in the 
figure label. Satellites with rg < 1 may have been affected by ''strangula- 
tion" . 



about subhaloes from the N-body, but track subhalo evolu- 
tion using semi-analytic recipes, we first compare the mass 
function and radial distribution of dark matter subhaloes 
predicted by the SAMs with the results directly obtained 
from the N-body simulations. We find that the parameteri- 
zations of subhalo merging and tidal stripping and destruc- 
tion adopted by the S08 and K08 models are able to fairly 
accurately reproduce simulation results for the mass distri- 
bution and radial distribution of satellites. 

We then test the luminosity function of our simulated 
satellite population against the latest observational results 
for the MW satellite luminosity function. Our models are all 
able to reproduce the LF down to a magnitude My = —5; at 
fainter magnitudes (—5 < Mv < —1) the K08 and S08 mod- 
els also provide a good fit to the observational data, while 
the MORGANA model tends to underestimate the abundance 
of ultra- faint satellites, though the predictions are still con- 
sistent with the MW data at the 1-a level. All models seem 
to suggest a decrease in the satellite number density below 
Mv ~ —5, consistent with the assumption of a NFW like 
radial distribution for observed satellites. 

We also perform the comparison between our model pre- 
dictions and the observations in the "observational plane", 
i.e., by applying "visibility" criteria to the simulated satel- 
lites and comparing with the SDSS data without any com- 
pleteness corrections applied. In this case, instead of assum- 
ing a radial distribution for the observed satellites (which 
could in principle depend on e.g. satellite mass or luminos- 
ity in a complex way) , we make use of the predictions of our 
models for the joint distribution function of satellite lumi- 
nosity and distance from the central galaxy. We again find 
good agreement, increasing the robustness of our results. 



We investigated the main physical processes responsi- 
ble for shaping the luminosity function of Milky Way satel- 
lites in our models. In the absence of cooling by molecu- 
lar hydrogen, and in the absence of processes like photo- 
ionization, SN feedback, or stellar stripping, the predicted 
satellite LF would show a peak at around Mv ~ — 14 and a 
sharp drop-off at Mv > —10, with essentially no satellites 
fainter than Mv ~ — 8 predicted. This drop below Mv > — 8 
is not due to limited numerical resolution, but rather to 
the sharp assumed cooling cutoff at T < 10 4 K (because of 
our adopted atomic cooling function). However, at temper- 
atures just slightly above the atomic cooling cutoff, cooling 
becomes quite rapid and so in the absence of some kind 
of feedback or suppression mechanism, these haloes rapidly 
cool all of the available baryons and convert them into stars. 

In our models, photo-ionization due to a cosmic reion- 
izing background and supernova feedback work together in 
order to re-shape this highly peaked luminosity function into 
the near-power law down to M v ~ —3 that is implied by the 
recent SDSS observations. Photo-ionization suppresses the 
infall of hot gas into low-mass dark matter haloes, reduc- 
ing the supply of baryons that are available for cooling and 
star formation, while SN feedback reheats cold gas and ex- 
pels it from small haloes, again suppressing the efficiency of 
star formation. If we include only photo-ionization or only 
SN feedback, we find an excess of intermediate luminosity 
satellites and a shortage of ultra-faint satellites (see also Ko- 
posov et al. 2009). In agreement with previous works (e.g. 
Kravtsov et al. 2004), we find that the satellite luminosity 
function in our "standard" models depends only weakly on 
the assumed redshift of reionization. 

We made use of results from the numerical hydrody- 
namic simulations of Gnedin (2000) to motivate our treat- 
ment of the suppression of gas infall due to the presence of a 
photo-ionizing background. A key parameter in this recipe 
is the "filtering mass", or the halo mass below which the 
gas content is significantly reduced relative to the cosmic 
average. However, recent work by other groups (Okamoto 
et al. 2008; Hoeft et al. 2006) has found that the filtering 
mass may be considerably smaller than the results of Gnedin 
(2000) suggested. We investigate the implications of modi- 
fying the normalization and/or redshift dependence of the 
filtering mass as suggested by these works, and find that 
when a lower normalization of the filtering mass is adopted, 
the results are more sensitive to the redshift of reionization. 
We find that we can still reproduce the observed luminosity 
function with the lower filtering mass if we adopt a higher 
reionization redshift (z r ~ 11 instead of z r ~ 7.5). 

We investigate the impact of stellar stripping on the 
observed luminosity function using the MORGANA model. 
We find that stellar stripping can only decrease the satellite 
stellar masses by at most about ~ 20%, and therefore prob- 
ably has only a minor effect on the satellite LF. However, 
we find that the modelling of tidal destruction of satellites 
does have a significant effect on faint end of the predicted 
LF. 

In this work, we have concentrated on the comparison 
with observations of one particular quantity, the statisti- 
cal distributions of satellite luminosity, although naturally 
the SAMs provide predictions of many other galaxy prop- 
erties (e.g. Lagos et al. 2009). This has been done mainly 
for two reasons: first because robust state of the art obser- 
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vational data have recently been made available for those 
two quantities and second because our goal was to directly 
address the so called missing satellite problem and the ori- 
gin of the newly discovered ultra-faint population. Other 
observed properties (such as metallicity, gas content, star 
formation history, etc.) contain complementary information 
on the formation mechanism of galactic satellites and more 
can be learned by studying in them in details (e.g. Okamoto 
et al. 2009). While the objects produced in our SAMs do 
resemble several properties of observed satellites we have 
decided to defer a more extensive explorations of those prop- 
erties to a future work. 

The semi-analytic models that we use in this work were 
originally normalized to reproduce global quantities such as 
the field galaxy luminosity function or gas fractions for rela- 
tively luminous galaxies (My < —16). They have not previ- 
ously been extensively tested against observations of galax- 
ies on these very small mass scales. In the case of the K08 
and S08 models, we found that the identical model ingredi- 
ents and even parameter values used in the standard versions 
of these SAMs (e.g. Somerville et al. 2008, Kang et al. 2005) 
were also able to reproduce the MW satellite LF. In the 
case of the MORGANA model, we found that the SN feed- 
back recipe in the original model (Monaco et al. 2007) had 
to be modified in order to reproduce the faint satellite pop- 
ulation. In either case, we have gained important insights 
about the physical recipes incorporated in these models and 
their applicability over a wide range of galaxy mass scales. 

Our study confirms and expands on previous works that 
address the so-called missing satellite problem using semi- 
analytic models and simulations. The main new contribu- 
tion of our paper is the implementation of SAMs within 
merger trees extracted from numerical simulations with very 
high mass resolution (particle mass ra v ~4x 10 5 fe _1 Mq), 
in which each of our four Galaxy-sized halo simulations 
contains 2-4 million particles. Thus, unlike previous stud- 
ies (e.g. Benson et al. 2002; Somerville 2002; Kravtsov 
et al. 2004) which compared only with the "classical" satel- 
lite population My = —9, we can resolve the very small 
subhaloes that may host the newly detected population of 
ultra-faint satellites. We also have the advantage of multiple 
simulations, unlike studies based on the Via Lactea simula- 
tion (Diemand et al. 2007) alone. Although our results are in 
qualitative agreement with studies based on simpler analytic 
recipes for assigning baryons to dark matter (sub)-haloes 
(e.g. Koposov et al. 2009), we showed that in our models, 
the MW satellite LF is shaped by a complex combination 
of different physical processes including tidal destruction, 
photo-ionization, and supernova feedback. 

Our models neglect several other physical processes that 
have been discussed in the literature, and which may be im- 
portant in shaping the properties of galaxies on these mass 
scales. Although we model the suppression of gas infall by a 
uniform cosmic radiation field after reionization, we do not 
account for the modification of the atomic cooling function 
by the radiation field, the possible photo-evaporation of gas 
from small haloes after reionization (Barkana & Loeb 1999), 
or photo-ionization by the nearby large galaxy (i.e., the ra- 
diation field of the Milky Way; e.g. Weinmann et al. 2007). 
Perhaps most importantly, we neglect cooling via molecular 
hydrogen, and the associated complex and poorly under- 
stood possible positive and negative feedback effects con- 



nected with the formation and destruction of H2 (e.g. Sal- 
vadori & Ferrara 2009, Ricotti et al. 2008). The fact that 
our models are nevertheless able to reproduce the bulk of 
the ultra-faint satellite population may indicate that these 
other processes operate at second order or cancel each other 
out, or may be simply a fortuitous coincidence. Certainly 
this bears further study. 

A further concern is that standard SAMs are known to 
fail to reproduce the more detailed properties of satellites 
in larger mass hosts: although they correctly reproduce the 
number of satellites as a function of halo mass, several dif- 
ferent SAM codes (including the ones considered here) have 
been shown to produce too large a fraction of red and pas- 
sive satellites (e.g. Kimm et al. 2008 and references therein) 
compared with observations. The main cause of this diffi- 
culty is believed to be the standard assumption that the hot 
gas halo is immediately stripped when a satellite enters a 
larger host, thereby depriving satellites of any new supply 
of gas (see Kang & van den Bosch 2008 for a detailed dis- 
cussion and possible solution). It is unclear whether this will 
impact our predictions for the very low-mass MW satellites 
— we defer this question to a future investigation. 

If we accept this "baryonic" solution of the "missing 
satellite problem" , other interesting implications follow. Our 
analysis predicts that roughly 1/5 of the total number of 
subhaloes with a present-day bound mass M > 2 x 1O 7 M 
should be dark. In order to properly test this picture, a 
signature of the presence of these dark satellites is needed. 
One possibility is that they could be detected via gravita- 
tional lensing (e.g. Metcalf & Zhao 2002) since those small 
subhaloes will act as perturbers of the lensing signal com- 
ing from the main halo. Unfortunately recent results based 
on numerical simulations have shown that perturbations in 
the lensing potential induced by (dark) satellites are very 
small and unlikely to explain the anomalous flux ratios of 
some multiple lensed QSOs (Maccio et al. 2006; Maccid & 
Miranda 2006). Another possibility would be detection of 
7-rays from dark matter annihilation, as the presence of 
substructure boosts the 7-ray signal by a factor of 4 to 15 
relative to smooth galactic models (Diemand et al. 2008). 
Finally, a third possible opportunity to detect the presence 
of a significant dark population of subhaloes in the MW halo 
could come from the signatures of the interaction of such a 
population with the thin stellar streams in the MW halo 
( e.g. Odenkirchen et al. 2001, Ibata et al. 2002, Johnston 
et al. 2002, Grillmair & Dionatos 2006). 

In final summary, our results show that not only is there 
no longer a "missing satellite problem" , but that well-known 
and well-motivated astrophysical processes working within 
the ACDM framework naturally predict the form of the ob- 
served MW satellite luminosity function over six orders of 
magnitude in luminosity. Indeed, it may be that convincing 
proof of the existence of the large predicted population of 
dark subhaloes via one of the methods suggested above (or 
one not yet discovered) is one of the last remaining major 
challenges for the ACDM paradigm. 
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